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ABSTRACT 

We examine the formation and evolution of the density enhancement (density 
spike) that appears downstream of strong, cosmic-ray-modified shocks. This 
feature results from temporary overcompression of the flow by the combined 
cosmic-ray shock precursor/gas subshock. Formation of the density spike is 
expected whenever shock modification by cosmic-ray pressure increases strongly. 
That occurence may be anticipated for newly generated strong shocks or for 
cosmic-ray-modified shocks encountering a region of higher external density, for 
example. 

The predicted mass density within the spike increases with the shock 
Mach number and with shocks more dominated by cosmic-ray pressure. For 
very strong shocks, the total compression compared to the upstream gas may 
approach ^^-D during the formation period, where 7^ is the gas adiabatic index 
and D is the compression ratio through the precursor. As the full shock reaches 
equilibrium, the spike detaches, lags behind the modified shock transition and 
is further compressed, so that the density can exceed the limit quoted above. 
We find this spike to be linearly unstable under a modified Rayleigh- Taylor 
instability criterion at the early stage of its formation. Our linear analysis 
shows that the fiow is unstable when the gradients of total pressure (gas 
pressure -|- cosmic-ray pressure) and gas density have opposite signs. We 
confirm this numerically using two independent codes based on the two-fiuid 
model for cosmic-ray transport. These two-dimensional simulations show that 
the instability grows impulsively at early stages and then slows down as the 
gradients of total pressure and gas density decrease. Flow within the density 
spike becomes disordered through the instability. It seems likely that this can 
significantly increase the local magnetic field beyond compressional effects. 
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Observational discovery of this unstable density spike behind shocks, possibly 
through radio emission enhanced by the amplified magnetic fields would provide 
evidence for the existence of strongly cosmic-ray modified shock structures. 

Subject headings: cosmic-rays — hydrodynamics — shock waves 
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1. Introduction 

Energetic charged particles (Cosmic-rays or CRs) are believed to be accelerated in 
strong shock waves through the first-order Fermi mechanism when they become trapped 
there by Alfven wave turbulence. A relatively small population of particles can achieve 
very high energies by repeatedly crossing the shock and having their velocities redistributed 
by scattering. This very efficient, diffusive acceleration mechanism was discovered by a 
number of people independently in the late 1970s ( [Krymsky 1977| , |Axford et al., 1977] , [Bell 



1978| , [Blandford fc Ustriker 19751 ). The process is nonlinear, since pressure feedback from 



accelerated CR particles modifies the shock compression and adjusts the acceleration 
efficiency towards an equilibrium {e.g., [Blandford 1980| ). The modified shock structure 



includes an extended upstream pressure gradient ("foot" or "precursor") that results from 
upstream diffusion of CR pressure. Within the precursor gas is adiabatically compressed. 
At least during formation of the precursor there will also be a dissipative, gas subshock. 
However, above some critical Mach number diffusive shock theory suggests that this feature 
can become weak or absent (leading to a smooth or "C-shock" structure) in CR shocks 
as they approach dynamical equilibrium {e.g., Urury fc I'alle l9^ ; [Jones fc Kang 199U ). 



The presence of a strong magnetic field transverse to the shock normal tends to reduce the 
efficiency of CR acceleration in shocks of a given sonic Mach number ( [Webb et al. 1986[ , [Jun 



et al., 19"9^ [Frank et al. 1995| ), so the critical Mach number for the formation of a C-shock 



will generally be higher for MHD shocks ([Jun et al., 1994[) . 



The hydromagnetic structures of CR shocks are now fairly well studied for a variety 



of circumstances. Much of that effort has focussed on steady-state conditions {e.g., [Drury 



fc Volk 198lt IWebb et ffl/.^986t [Ellison, Mobius fc Paschmann 199d|; [Kang fc Jones 



1990 ; Jun et al., 1994[) . On the other hand, there are important circumstances in which 



steady-state CR shocks may not be applicable, such as supernova blast waves {e.g., porfi 
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1984; Drury, Markewiecz fc Volk 198S ; Jones fc Kang 1990 ). One interesting feature in 



time-dependent, evolving CR shocks is the appearance of a sharp density enhancement or 
"spike" downstream of the shock transition in both parallel and perpendicular shocks (Porfi 



|1984| ; Porfi 1985| ; [Jones fc Kang 1990| ; |Jun et al., 1994| ). This sharp density enhancement is 



formed as the growth of CR pressure first enhances the total compression through the shock 
and then levels off towards the equilibrium predicted for steady-state shocks. The maximum 
density within the spike can exceed that for the equilibrium modified shock and the initial 
shock by a large factor when the shock Mach number is large. We have discovered that this 
feature is Ray leigh- Taylor unstable. That could have interesting observable consequences, 
if, for example, it lead to amplification of magnetic fields within a well-defined region 
behind a supernova shock. In this paper, we present one- and two-dimensional numerical 
simulations of the evolution of the density spike, demonstrating its instability. We point out 
that the instability discussed here is distinct from other previously reported instabilities 
associated with the CR shock structure itself. On the other hand, the instability requires a 
large CR pressure and occurs because of modifications to the shock. Our simulations are 
based on the two-fluid formulation for CR transport ( prury fc Volk 1981| ). This method is 



much simpler than the more complete diffusion-convection formalism of CR transport, and 
it is still the only really practical method for studying multidimensional, time- dependent 
CR modified flows. Although the two- fluid method has limits {e.g., [Jones fc Ellison 1991|) , 



recent work has demonstrated that the two-fluid formalism gives consistent solutions 
with diffusion-convection formalism in time dependent shocks when the necessary closure 
parameters are carefully deflned [e.g., [Kang et al., 199^ ; prank et al. 1995t [Kang fc Jones 



1995t [Kang fc Jones 1996| ). 



The outline of this paper is as follows. The detailed description of the density spike 
formation is given in §2. In §3, we study the formation and evolution of the density 
spike in a various shock structures by means of one- dimensional numerical simulations. 
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The modified criterion of Rayleigli- Taylor instability in the presence of CRs is discussed 
by solving the linear perturbation equations of the two-fluid system in §4. §5 presents 
our two-dimensional numerical results of the unstable density spike simulation, and we 
summarize in 56. 



2. Formation of the Density Spike 

Formation of the sharp density enhancement in CR modified shocks has been 
explained previously ( Prury 1987| , [Jones fc Kang 1990| ). We describe the process of density 
enhancement within CR-modified shocks more fully here in order to understand better 
what conditions may be relevant to the new instability we report. The sharp density 
enhancement or "spike" originates from development of or increases in the CR precursor. 
Diffusive acceleration at shocks depends on the ability of a fraction of high energy particles 
isotropized by scattering in the postshock flow to return upstream of the shock. Most of 
those particles are then advected back through the shock in a process that repeats. On 
balance these two effects confine the upstream CRs within a region whose scale is the 
diffusion length, ~ i^/vq, where k is a suitably defined mean spatial diffusion coefficient 
Drury fc Volk 19MD , and Vq is the upstream flow speed of the gas relative to the shock. 



see 



If we define Xd in terms of the cosmic ray pressure, Pc, and the CR pressure increase before 
the shock is 6Pc, then there is a pressure gradient in front of the shock, ~ 6Pc/xd- That, in 
turn, decelerates the gas flowing into the shock, which produces an adiabatic compression 
and heating of the gas, adding to the total pressure gradient. 

For strong, high Mach number and steady shocks this "precompression" and 
"preheating" can be sufficient to eliminate the need for a dissipative, gas subshock in order 
to satisfy the steady jump conditions for the shock {e.g., [Drury fc Volk 1981| ; [Malkov fc 



^olk 1996[). Then, according to this theory, the steady transition between the upstream 
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and downstream states will be continuous (a C-shock is produced). On the other hand 
such a shock will usually form from one that initially is not dominated by CR pressure. In 
other circumstances the Mach number of the shock may increase, for example, as the shock 
encounters a denser, colder environment or as the shock accelerates down a steep density 
gradient. That shock strengthening is followed by an analogous increase in the CR pressure 
and readjustment of the shock structure. In either of these two situations even an eventual 
C-shock will pass through an intermediate condition with a precursor and a gas subshock. 
In this transitional stage, the fluid is compressed twice; once by the precursor and once by 
the gas subshock. The precompression of the gas within the precursor automatically reduces 
the Mach number of the subshock. However, if the subshock is initially strong, then while 
the CR pressure at the shock is small compared to PqVq, the second density compression, 
through the subshock, is hardly reduced from its initial value. Consequently, the double 
transition of the flow can result in a density compression of the flow much more than the 
usual times the density far upstream, where 7^ is the adiabatic index of the gas. 

In order to determine approximately the maximum possible compression ratio without 
solving the full set of equations we can just use the above double compression concept. Let 
us define four distinct states in the flow for convenience. State "0" refers to flow conditions 
far upstream, ahead of the precursor. State "1" is that just upstream of the gas subshock, 
if it exists, while state "2" is that just downstream of the subshock. Then, finally "3" refers 
to the state that existed downstream of the subshock before the shock began adjusting to 
new upstream conditions. That region will propagate downstream, away from the subshock, 
and it is in the space between states "2" and "3" that the density spike feature develops. 
We also define the compression ratio through the precursor as D = ^. Recall that in this 
smooth transition the gas will be compressed adiabatically, so that Pg oc p^" . Then, one can 
use the gas adiabat and, assuming that the compression through the precursor is slowly 
changing, mass flux conservation to obtain the relation between the sonic Mach numbers 
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(M) in states "0" and "1". 



Ayf _ ^1 _ ^0 ^ _ Mq , 

■'^1 ~ ~ 73 + 1 ~ 7g + l 

ci Co 



where c is the gas sound speed. Therefore, we can see that the sonic Mach number of 
the flow decreases as the flow experiences the adiabatic compression. If Pcs is less than 
the equilibrium value, then P^i = Pco + ^Pc will increase, causing D to rise as well. The 
value of D at equilibrium in a given shock depends on various details such as the Mach 
number of the shock, but, from conservation laws cannot exceed D = D^ax = 4 if the net 
adiabatic index of the combined thermal/CR plasma is nonrelativistic or D = Dmax = 7 if 
that plasma is relativistic. These limits would correspond to strong shocks that have been 
smoothed entirely into C-shocks by CR pressure ( prury fc Volk 1981| ). 



As D increases for a given M.q the value of A^i will decrease, down to a limit M.i = 1. 
If the entire steady-state jump conditions can be satisfied when A^2 > 1 the shock should 
evolve into a C-shock in which the downstream fiow is supersonic ( prury fc Volk 198 1| ). If 
the steady jump conditions across the full transition require a subsonic fiow in region "2", 
then a subshock is required even as the shock reaches equilibrium. Applying the condition 
M.I = M.2 > 1 for D = Dmax leads to a minimum constraint on the Mach number Aio for 
a CR shock to be able to evolve into a C-shock; namely 

7g + l 

Mo > Dnix ■ (2) 

The critical Mach number for a C-shock lies in the range 6.35 to 13.39 for Dmax between 
4 and 7; that is as the net adiabatic index changes from 7c = 5/3 to 7c = 4/3. It should 
be noted that this critical sonic Mach number is approximate because D in C-shocks is in 
general slightly smaller than Dmax- Therefore, the critical sonic Mach number for C-shock 
should be slightly smaller than our estimated value. The critical sonic Mach number 
increases for perpendicular shocks when the magnetic field is strong (|Jun et al., 1994| ). In the 
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perpendicular shock case, the condition for a C-shock is fi(= f|) > c^(= c^) + v\ ^^2) 
where is the Alfven speed. This gives the critical sonic Mach number 



^2 1/2 

•Mo > + 4^^LJ • (3) 



As long as the gas subshock remains, we can simply relate the states "1" and "2", by 
the standard shock jump conditions. An assumption in diffusive shock theory is that there 
is no discontinuity in Pc, since the CRs necessarily are able to cross the subshock freely. If 
there is a subshock, we can write 

pi (7, - l)M\ + 2 ( (Igii)^ + 2)D7.+i • ^ ^ 

If the subshock has vanished p2 = Pi- Including a subshock the total density jump between 
states "0" and "2" is then 

P2 ^ {% + l)Ml 

If we assume that the Mach number A^o is large compared to the critical Mach number in 
equation 2, but there is still a strong subshock equation 5 reduces to 

Ei^JiHd (6) 

Po 79-1 

This last expression cannot apply to a steady shock, but can be approximately true during 
the evolution of the shock towards its time asymptotic state. For a relatively brief time 
compared to the evolution time scale of the shock it can even happen that the subshock is 
strong (so that A^i 00 in equation 4 but that D ^ Dmax)- Then the total compression 
can exceed that of either the precursor or the subshock separately. In fact it can lead 
to compression as high as 28 (16) for 7^ = 4/3 (7c = 5/3) just as the shock approaches 
an equilibrium and the spike detaches (the result indicated by equation 6). The extra 
compression is also accompanied by a modestly enhanced total velocity jump during this 
time {e.g., [Jones fc Kang 1992| ). This, in turn, can increase the rate of particle acceleration. 
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at least for high energy particles whose scattering lengths are large enough to allow them to 
cross substantial parts of the precursor each time they pass through the subshock ( Prury 



1987|) . This influence can speed up the modification process itself as seen in Fig. 5 of |Jun et 



al., 1994 



The overcompression forms a spike, because pa, if it represents an equilibrium condition 
for the earlier shock, will satisfy the constraint, = ps/po < Dmax- Hence, for the reasons 
we just discussed, there is a time when p2 = r2Po > Ps- We can estimate the width of the 
spike through a simple consideration of that transition interval. If the time to reestablish 
equilibrium is teq, then the width of the spike is Xg = V2teq, where V2 is the mean velocity 
relative to the shock inside the density spike. We can assume that the density p2 peaks in 
time just before the shock reaches its steady structure and define the total compression at 
that time to be = p2{peak) / pq. Simplifying the density structure of the spike to a linear 
profile we can estimate V2 ~ r3+r2 ' '^^^^ '^^^ ^ simple analytic estimate for t^q given 
by Jones & Kang (1990) in terms of the upstream diffusion time, td = Xd/uo = k/u^, the 
ratio, Pc2/Pco of the equilibrium downstream and upstream cosmic-ray pressures and jc] 
namely, 

5td, MPc2/Pc0) + l 



te.^yln r---7 1 , (7) 



where g = — 1. In the limit (7 — > 0, teq ~ 5td{Pc2/ Pco — 1)- This leads to a simple 
estimate for the width of the density spike at "separation" ; namely, 

9{r3 + r2p) g + l 

with an appropriate limit taken as (7 ^ 0. For moderate to strong shocks this formula will 
generally predict Xs/xd ~ few x(10 — 100), with larger values for larger Pc2/Pco and as 
7c ^ 4/3. 
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3. One-Dimensional Numerical Simulations 

In order to study the formation process and evolution of the density spike in 
CR-modified shocks, we use the modified-ZEUS code which solves two-fluid equations on 
an Eulerian mesh. The algorithms of the code are described in Jun, Clarke, and Norman 
(1994). 

Figure 1 shows the evolution of CR-modified shocks and density spikes computed 
for different total Mach numbers defined by vq/ \J'^^^ + -^^^ To simplify our discussion 
we have chosen a constant diffusion coefficient, k = 1 for the moment. The value of the 
diffusion coefficient enters primarily through its role in determining the width of the density 
spike (equation 8) and through its role in controlling the rate at which the CR modifications 
to the shock take place. Figs, la, lb, and Ic show the evolution of CR-modified shocks with 
Mach number 10, 20, and 50 respectively. The adiabatic index for CRs is 7c = 5/3, for these 
particular tests. Other values will be included later. A shock wave is initially generated by 
supersonic fluid striking a reflecting boundary at x = and moves to the right. The dotted 
lines represent the shock structure at the earliest stage while the later stages are represented 
by the dashed line and solid lines in a time sequence. The density overshoot is clearly seen 
in early stages of evolution. This density overshoot results in a "density spike" as the shock 
evolves further and the gas subshock becomes weaker due to the decreased gas pressure 
which is replaced by CR pressure. The shock front is strongly modified by the CRs and the 
gas subshock disappears at later stage of the evolution so that the shock transition becomes 
continuous (C-shock). The formation of the density spike occurs as the density overshoot 
lags behind the shock. This lagging starts after temporary high pressure (overshoot) from 
the shock front generates a positive total pressure gradient and starts to push the density 
spike downstream. As the overcompressed density spike moves downstream, the shock 
structure reaches a steady state that can be determined by conservation equations ( prury fc 
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Volk 1981, Webb 19'83 , Jun et al., 1994| ). The density spike is a moving compression region, 



and material flows through it. We point out that this feature is not a contact discontinuity. 
The total compression ratio of the density spike is lower than 16 (for 7c = 5/3) when it 
is formed at the shock front. However, after the density spike moves downstream, the 
compression continues and increases the density further so that it can become higher than 
16 for 7g = 7c = 5/3 (see Fig.lc). The flow compression eventually slows down due to the 
finite pressure gradient and increasing density. We find in the overcompressed region that 
the density becomes higher in a stronger shock (higher Mach number). The development 
of the density spike is not restricted to the piston driven shock problem. Figs. Id, le, 
and If show the evolution of CR-modified shocks which were initialized as a standing 
shock by assigning analytic values for upstream and downstream states. Total pressure in 
the downstream is divided equally to gas pressure and CR pressure. The left boundary 
condition is outflow and the right boundary is inflow with an upstream condition. The 
results show qualitatively the same density spike evolution as for the CR-modified piston 
driven shocks. Once detached, the density spike appears to move leftward since the shock 
front is at rest. The density spike becomes sharper as Mach number increases because the 
diffusion length is inversely proportional to the shock speed (see equation 8). 

In Fig. 2, we present the evolution of three different CR-modified shocks. Fig. 2a shows 
the formation of the density spike in a Mach 20 standing shock with the CR adiabatic index 
7c = 4/3. The shock transition becomes continuous eventually due to the CR domination 
in the downstream state of the shock. This CR domination of the shock also results in the 
steady-state density jump of the shock higher than 4. The prominent density spike is also 
found in the downstream of the shock as expected. Fig. 2b shows the result of CR-modified 
shock with the diffusion coefficient k, = 1/p. The rest of initial conditions are identical to 
the case presented in Fig. Id. As reported by Jun et al. (1994), the steady state is obtained 
faster than in the case with k, = 1. Another feature is the sharper structure of the shock 
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and density spike because of the smaller diffusion length in the shock transition region 
and downstream state where the density is higher than 1. Fig. 2c shows the evolution of 
CR-modified shock which contains a gas subshock and a smooth precursor. The density 
spike is also clearly present, although it is not as prominent as in C-shock case. This 
example manifests that the density spike is not restricted to the C-shock case. This density 
enhancement always follows the formation or enhancement of a precursor. Therefore, 
the density augmentation should form whenever the CR acceleration is efficient and the 
back-reaction of CR to the shock front is dynamically important. 

The density spike also forms in simulations involving shock propagation into nonuniform 
media. It has been seen in cases where the shock moves down a steep density gradient, such 
as a SNR blast into a preexisting wind cavity, or where a CR-modified shock encounters 
a steep increase in upstream density, such as a cloud. Fig. 3 shows the evolution of a 
Mach 10 shock propagating into a higher density background. The shock is initialized as 
a standing shock at x=3. The shock is already modified by generating a precursor and 
density overshoot at t=0.1. By t=0.2, a density spike has formed and is moving to the left 
relative to the shock front. The shock transition has reached a steady state and the gas 
subshock disappeared (leaving a C-shock). At t=0.3, the C-shock is approaching a higher 
density medium (p = 3), while at t=0.4, the C-shock is interacting with the higher density 
region and another density overshoot is produced. Finally, at t=0.5, the shock reached to 
the new steady state and the second density spike is formed. The C-shock structure is 
recovered (Figs. 3e and 3f). 

As the CR-modified shock encounters the denser external medium, its precursor is 
compressed and the shock structure temporarily reforms a viscous, gas subshock. This 
process can occur if the rate of increase of the incoming gas density exceeds the rate 
at which the shock is able to adjust to an equilibrium in terms of its CR pressure, as 
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given by equation 7. As the newly formed discontinuous shock propagates into the denser 
background, the shock is modified anew and the precursor is enhanced. The important 
efi^ect during the shock passage into the denser region is a higher acceleration efficiency due 
to increased compression. If the inhanced upstream density is in pressure equilibrium with 
its surroundings, it will be colder. Then it is easy to show that the Mach number of the 
shock increases as it penetrates the denser medium. The importance of this will be more 
significant in an already CR-dominated C-shock, because the full shock compression can 
become much higher in responce to the density enhancement encounter than through the 
original incoming C-shock. Jones and Kang (1992) found in their shock-cloud two-fiuid 
simulations that the acceleration efficiency is increased more by the encounter if the incident 
shock is previously dominated by CR pressure. 

In the begining of the density spike formation, the total pressure decreases upstream of 
the density spike in order to conserve total momentum fiux {pv^ -\- Pg-\- P^. As a result, the 
total pressure in the left side of the density spike has a negative gradient while the density 
gradient is positive. As we will show in the next section, these opposite gradients between 
the total pressure and the density consitute an unstable condition. As the shock structure 
itself reaches a dynamical steady state (equilibrium), the momentum fiux upstream of 
the density spike becomes constant and the steepening of the density spike slows down. 
Eventually, the gradients of the total pressure and the density decreases, after the density 
spike is detached from the shock front, mainly because of CR diffusion. Once this happens, 
the growth rate of the instability decreases accordingly. Prom this dynamical picture it 
is obvious that the duration of the instability is comparable to the shock equilibrium 
timescale (equation 7) each time the shock readjusts. Since the opposite gradients of the 
total pressure and the density exist only for a short period of time during the formation of 
the density spike, the driving force of the instability appears to be impulsive. 
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4. Modified Rayleigh- Taylor Instability 

One can expect that cosmic-ray pressure should be added to the total pressure in the 
criterion for the Rayleigh- Taylor instability in the presence of cosmic-ray particles. In this 
section, we derive the modified criterion for the Rayleigh- Taylor instability in the presence 
of CR pressure by considering the two- fluid model (Prury fc Volk 1981|) . The governing 
equations are: 

f^+V-{pv) = (9) 

dv 

p— + piv ■ W)v + WiPg + P,)+pg = (10) 
dP 

-^+v.VP, + ^,P,V-v = (11) 

dP 

—^ + v.VP, + 7cPcV ■v = V- (kVP,), (12) 
ot 

where Pc is the CR pressure and k, is the mean spatial CR diffusion coefficient. The gravity 
term is added to the momentum equation to satisfy the equilibrium state in a stationary 
condition. For this analytical analysis we assume k to be a constant, but that is not 
necessary, either for this instability, or for validity of the two-fluid model (see [Kang et"ar 



19921 ; [Kang fc Jones 1995| ; [Kang fc Jones 1996] for specific tests of time-dependent two-fluid 



simulations against more detailed diffusion-convection equation or "particle" simulations). 
We consider an exponential background, 

% = Voexp{ ^ ^° ), p = poexp C j Pg = Pgoexp{^^——^), Pc = Pcoexp C ^ ^° )- 

Li L2 L3 L4 

(13) 

In general, the normal mode for a nonuniform background contains a dependence on the 
x-direction, and one needs to solve the characteristic value problem with proper boundary 
conditions. However, since we are seeking the instability criterion only, we can restrict 
our analysis to the local region, x — Xq <^ Li. We consider plane wave perturbations, 
~ ipiexpint + ikyn) where ipi is a perturbed quantity {v^i, Vyi, pi, Pgi, or Pd). Then, one 
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can obtain the linearized perturbation equations using equations 9-13, 

piin + ^) + ikypoVyi + y^t'xi = 0, (14) 

v^i{n^—)--^{—^—) = ^, (15) 

ni;,i + ^(P,i + Pei) =0, (16) 
Po 

^91 (^ + 73 -p-) + t-xi^ + iky-igPgoVyi = 0, (17) 

Pciin + 7c + k/cJ) + fxiT^ + iky'jcPcoVyi = 0. (18) 
By combining the above equations, we get the dispersion relation, 

X {n{n + ^)(- + ^ + -kl) + kla%n + ^ + nk^ + ky^ + ^)) 

- S^(n + ^ + kA:J)(^ + ^) - S^(n + ^)(^ + ^) = (19) 

PO-^3 -^1 ■ -^3 -^4 Po-^4 ^1 ^3 J^A 

where a, = (l£^)i/2 and = (2£^)i/2. 

Let us consider two different limits of this dispersion relation in a stationary background 
{vo ^ 0). 



4.1. Short-Wavelength Limit (ky — > 00) 

The familiar and related Ray leigh- Taylor and convective instabilities in pure gas 
dynamics are actually somewhat different in terms of criteria, because the convective 
instability depends on fluid compressibility, whereas the classical Ray leigh- Taylor instability 
does not. A flow is convectively unstable if ^^/^^^ < l/Tg? while the flow is Rayleigh- 
Taylor unstable if §^/^l:2Es. < g (Bandiera 1984). This shows that the convective instability 



is present in a wider range of conditions than the Rayleigh- Taylor instability. When the CR 
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diffusion term is ignored, our dispersion relation reduces to the growth rate for the general 
convective instability in the short-wavelength limit (for the frequency range, n ^ nky but 

pl(al + al) ^Ls^ lJ p^L^ ^ Ls ^ lJ ' ^tPo^ dx dx ^ 

where 7* = ^^^^o+Z-^-" and Pt = P„ + Pc- 

In the presence of CR diffusion and short wavelength limit (for the frequency range, 
n <C nky and <^ ky{ag + a^), the dispersion relation reduces to 

= + ^)(Jl 'i. ]_dP^AdlnP^ _ dbwy .gl) 

po -L3 L4 ^gLs L2 po dx 7g dx dx 

This relation includes criteria for both gas convection and Rayleigh- Taylor instabilties. 
Now, let us consider our dispersion relation for several different cases. 

i) First, when all three gradients {dPt/dx,dPg/dx, and dp/dx) have same sign, then 
the flow can be convectively unstable if l^^^^l > 1^1. 

1 7g ox I ' ax ' 

ii) Second, the flow is always Rayleigh- Taylor unstable if 

f >0,f >0,|^<0. or f <0,^<0,|^>0. (22) 
dx dx dx dx dx dx 

iii) Third, the flow is Rayleigh- Taylor unstable if 

^>0,f <0,|2<0, or f <0,f >0,|^>0, (23) 
dx dx dx dx dx dx 

when the gas is convectively stable ( | ~ '^^q^^ I < I^^D- ^ is interesting to see that the 
convective instability can be suppressed by the opposite CR pressure gradient. 
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4.2. Large- Wavelength Limit {ky — > 0) 

In the large- wavelength limit (for the range, n ^ nky and ^ ^1(^1 + '^c))) 
dispersion relation reduces to 



= _^(^ + ^) ^ 1 gp ^(P. + Pc) _ .24) 
L2P0 -^3 -^^4 Po dx 

Therefore, the flow is unstable if 



dp dPt ^ ^ ^25) 
dx dx 

This is the generalized crossed-gradient condition for Rayleigh- Taylor instability when the 
total pressure is considered instead of the gas pressure (e.g. see [Jones et al., 1981| ). One can 
see that the convective instability disappears in the large- wavelength limit. 

From the above analysis, one can see that the flow is Rayleigh- Taylor unstable 
throughout the entire wavelength range when the total pressure and gas pressure gradients 
are in opposite sign to the density gradient, which is the instability condition in our density 
spike case. 



5. Two-Dimensional Numerical Results and Discussion 

The nonlinear development of the unstable density spike is best studied by means 
of multi-dimensional numerical simulations. We have simulated the instability of the 
density spike in two dimensions by setting up several different CR-modified shocks, such 
as piston-driven shocks and standing shocks. The density spike is found to be unstable 
in every case. To reduce the potential that we are observing only a numerical effect, we 
carried out analogous simulations with two very different codes; namely a version of ZEUS, 
modified to include two-fluid CR physics ( |Jun et al., 1994| ) and a version of PPM, also 
modified and used extensively in the past for a number of CR applications (e.g., [Jones 
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& Kang 1990; |Jones fc Kang 1993| ; [Ryu, Kang fc Jones 1993| ). The gasdynamical methods 



are quite distinct in the two codes, and the CR transport methods are also different (one 
being exphcit and one being imphcit, for example). The results of the two codes were 
qualitatively the same, and most importantly, both demonstrated the instability. Fig. 4 
shows the two-dimensional results of a Mach 10 CR-modified shock structures computed 
by the modified ZEUS code. The one-dimensional time-dependent evolution of the same 
shock is shown in Fig. Id. The computational space is resolved by 600x200 uniform zones. 
A random density perturbation (10%) is carried in by the flow from the right boundary. 
The density perturbation (pressure-free perturbation) is chosen because we do not intend 
to generate sound waves that are unstable in the CR-dominated shock (acoustic instability, 
see, e.g., |Drury fc Falle 1986| ; [Kang et al., 1992]) . In this way, we should be free of any 



influence from the propagating acoustic instability, so that the instability we observe in the 
density spike results only from the Rayleigh- Taylor instability. 

Images in Fig. 4 show the density at t=0.1, 0.2, and 0.3 (top to bottom). A number of 
thin fingers are produced in the left side of the density spike as a result of the instability. 
In addition, these fingers show unstable curly structures that are likely to be the result of 
the secondary Kelvin-Helmholtz instability which was triggered by the shear between the 
finger and the background flow downstream of the Rayleigh- Taylor unstable region. The 
evolution of the unstable flow is found to be different from the classical Rayleigh- Taylor 
instability in the following way. In the classical Rayleigh- Taylor instability, short wavelength 
modes appear first in the linear regime, because their growth rates are higher. However, 
the nonlinear Rayleigh- Taylor instability is generally dominated by larger wavelength 
features, because they reach higher terminal velocity and overtake smaller structures ( |Jun 



et al., 19"95| ). However, short wavelengths are still dominant in the nonlinear stage of our 



CR-modified shock case. This is because the unstable condition exists only temporarily and 
does not provide sufficient time for growth of perturbations on the scale of the thickness 
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of the density spike itself. (Recall that the perturbation wavevectors are transverse to the 
flow direction.) Short wavelength features are just stretched along the flow direction by the 
initial acceleration to form long thin fingers without much interaction with one another. 

To provide a sense of the time evolution of the instability, Fig. 5 shows the amplitude 
of the density fluctuations as a function of time. The amplitude, A(x) is deflned as 
jj Jo Sp{x,y)^dy where H is the width of the computational plane in the Y-direction and 
Sp{x,y) = ^ - 1. This amplitude is averaged over the region X2 — Xi in the unstable 
flow; that is, A = J.^^ A{x)dx. The domain for the average is chosen to be 10 zones 

leftward from the peak density in the density spike. The amplitude, A provides us some 
information about the overall density fluctuation history. We flnd that Fourier amplitude of 
the density perturbation is not a good way to measure the growth of the instability because 
the flow is moving through the unstable region. 

The apparent behavior of A(x) allows us to divide the growth history of the instabihty 
into 4 different stages as indicated in Fig. 5. In stage 1, the flow carrying density perturbation 
has not reached the density spike. The amplitude of the density perturbation increases 
suddenly in stage 2 due to the incoming density perturbation. The density perturbation 
grows while the flow with the density perturbation passes through the unstable region 
(stage 3) . Recall that the growth of the instability exists only temporarily while the flow 
in the left (downstream) side of the density spike is unstable. The, flnal, relaxation phase 
of the instabihty is shown in stage 4. The sudden decay of the instabihty growth is hkely 
due to the combined effects of relaxation of the unstable region and the nonlinearity of the 
instability. The instability is found to result in a slight reduction in density in the spike 
compared to the one-dimensional behavior at this time. However, this reduction is small 
because of the short duration of unstable flow at the early stage. 

Although our perturbations were designed to avoid generation of upstream sound 
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waves that might contribute to this instabihty, we recognized that passage of the density 
perturbations through the shock might also produce sonic perturbations. To avoid any 
possibihty that this might modify the instabihty of interest, we have done a simulation 
with the diffusion coefficient, n = 1/p, because it is known that CR-modified shocks are 
stable against the accoustic instability when k takes this form (Prury fc Falle 1986|) . Fig.6 



shows a greyscale image of a two-dimensional numerical result for the case with n = 1/p 
at t=0.2 (one- dimensional result for this case is presented in Fig. 2b). The density spike is 
still unstable, thus, confirming the existence of the physical Ray leigh- Taylor instability in 
the density spike. The shorter density fingers compared to Fig. 4b are due to the narrower 
density spike caused by the reduction in diffusion length and equilibrium time when n = 1/p 
(compare Fig. Id and Fig. 2b). In the presence of sonic perturbations , [Ryu, Kang fc Jones| 
|1993| discovered in two-dimensions there is a secondary, Rayleigh- Taylor class instability 
associated with the nonlinear accoustic instability, which generates accoustic turbulence 
within the shock precursor. Since that instability will also act as strong CR-modified shocks 
come to equilibrium, it may be natural to expect that the density spike will often be subject 
to strong perturbations. 

An obvious place where one can expect the unstable density spike to be formed may 
be in the supernova remnant blast wave structure. Numerical studies of CR-modified 
supernova remnant shocks reveal the appearance of the density spike in the evolution ( |Jones 



fc Kang 1992| ). Stronger modification of the shock structure by CR pressure was found 



in the simulations when the supernova remnant shock propagates into an inverse square 
density medium such as a wind-blown background. Therefore, this density spike could be 
found in the early stage of supernova remnant evolution while the supernova shock is still 
inside the presupernova wind. Discovery of this unstable density spike and consequent 
turbulent flow would provide strong, in situ evidence within supernova remnants for the 
action of diffusive shock acceleration and, especially for strong shock modiflcation by 
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energetic particles. We note finally that, although these effects are strongest for C-shocks, 
the development of C-shocks is not necessary for what we have discussed in this paper. It is 
necessary, however, that CR pressure become the dominant means of decelerating the flow 
through the shock. 

6. Summciry 

We have given a detailed description and numerical study of the density spike 
formation downstream of a CR-modified shock structure. The density spike originates 
from a temporary overcompression while the gas shock is modified by the CR pressure. 
The resulting compression ratio is found to increase with the flow Mach number. We 
emphasize that the density spike is a common relic of the shock modification by the CR 
pressure. The density spike is found to be Rayleigh- Taylor unstable at the early stages of 
its evolution. Our linear analysis shows that the flow is Rayleigh- Taylor unstable if the 
total pressure gradient and gas density gradient are of opposite signs. We have simulated 
this unstable density spike by solving the time dependent two-fluid equations numerically 
on a two-dimensional grid. It is found that the instability grows impulsively in the early 
stages, but slows down as the density spike begins to relax its structrure. The numerical 
results show a number of thin flngers growing downstream from the density spike. However, 
the mass density in the spike is decreased only slightly as a result of the weak instability. 

The simulations were performed on a Cray C90 at the Minnesota Supercomputer 
Center. The work reported here is supported in part by NSF grant AST-9318959 by NASA 
grants NAGW-2548 and NAG5-5055 and by the Minnesota Supercomputer Institute. 
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Fig. 1. — Density evolution of CR- modified shock structures for various Mach numbers. 
The shock is propagating into a uniform background of p = l,Pg = l,Pc = 1. Adiabatic 
index, 7, is assumed to be 5/3 for both gas and CR and the diffusion coefficient, k = 1. 
Each plots (dotted line, dashed line, and solid line) are represented in a time sequence. 
The computational domain is resolved by 300 uniform grid zones. Readers are referred to 
Dorfi(1984,1985), Jones & Kang(1990), and Jun et al. (1994) for the evolution of other 
dynamical variables, (a) Piston driven shock with Mach 10. The shock is propagating 
with the velocity Vg — 18.257. Three density plots are taken at t= 0.1,0.2, and 0.3. (b) 
Piston driven shock with Mach 20 {vg = 36.515). Three density plots are taken at t= 
0.033, 0.066, and 0.099. (c) Piston driven shock with Mach 50 (v, = 91.287). Three 
density plots are taken at t=0.008, 0.016, and 0.024. (d) Standing shock with Mach 10 
{vs = 18.257). The upstream (p = 1.0, t; = -18.257, Pg = 1.0, Pc = 1-0) and downstream 
(p = 3.884, V = —4.7, Pg = 124.75, Pc = 124.75) states of the flow are initialized with analytic 
jump condition of CR-modified shock. Total pressure in the downstream is divided equally 
to gas pressure and CR pressure. Three density plots are taken at t=0. 1,0.2, and 0.3. (e) 
Standing shock with Mach 20 (upstream state: p = 1.0, v = —36. 515, = 1.0, Pc = 1.0. 
downstream state: p = 3.97, v = — 9.197, P^ = 499.75, Pc = 499.75). Three density plots 
are taken at t=0.033, 0.066, and 0.099. (f) Standing shock with Mach 50 (upstream state: 
p = 1.0, Vs = -91.287, Pg = 1.0, Pc = 1.0. downstream state: p = 3.995, = -22.847, Pg = 
3124.7, Pc = 3124.7 ). Three density plots are taken at t=0.008, 0.016, and 0.024. 

Fig. 2. — (a) Standing CR-modified shock (Mach 20, Vg = 34.641). The adiabatic index 
for CR 7c = 4/3 and k = 1. Plots are taken at t=0.1 (dotted line), t=0.2(dashed line), and 
t=0.3(sohd fine). The numerical resolution of the entire computational space is 300 zones, 
(b) Standing CR-modified shock (Mach 10) with the diffusion coefficient k, — 1/p. Other 
physical conditions are identical to the case in Fig. Id. (c) Standing CR-modified shock 
(Mach 5, Vs — 8.66,7c = 4/3). The computational space is resolved by 800 grid zones. 
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Fig. 3. — Mach 10 shock {vg — 18.257) is propagating into the higher density region [p — 3). 
The adiabatic index for CR is 7c = 5/3, and k, — 1. The computational space is resolved by 
400 grid zones. 

Fig. 4. — Grey-scale images show the density of two-dimensional numerical result of CR- 
modified standing shock with Mach 10 at t=0.1, 0.2, and 0.3 (from top to bottom). The 
computational space is resolved by 600x200 uniform grid zones. One-dimensional result is 
shown in Fig. Id. 

Fig. 5. — The growth of the density perturbation in the density spike in Mach 10 standing 
shock (fig. 4). 



where the average is taken over the region X2 — xi in the unstable flow. 

Fig. 6. — Grey scale image of CR-modified shock generated by Mach 10 shock at t=0.2. 
The diffusion coefficient used is 1/p. Other physical parameters are identical to the ones 
used Fig. 2b case. The computational space is resolved by 600x200 uniform grid zones. 




where H is the width of the computational plane in Y-direction, 6p{x,y) 
p{x) is the averaged density over Y-direction. 



p{x) 



— 1, and 




(a) t=0.1 



(b) t=0.2 



(c) t=0.3 






higher 
density 
backgrourjd 



1.0 2.0 3.0 4.0 

X 

(d) t=0.4 



0.0 1.0 2.0 3.0 4.0 

X 

(e) t=0.5 



2.0 
0.0 



continuous 
shock 



0.0 1.0 2.0 3.0 4.0 

X 

(f) t=0.6 




20.0 
15.0 



I 10.0 

0) 
■D 



0.0 1.0 2.0 3.0 4.0 

X 



5.0 
0.0 



1 1 1 1 1 

density spike 


: \ 






\ shock 

y : 








,1 



0.0 1.0 2.0 3.0 4.0 

X 




0.0 1.0 2.0 3.0 4.0 

X 



